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Abstract 

Predicting the dynamic behavior of a large network from that of the composing modules is a central problem in systems 
and synthetic biology. Yet, this predictive ability is still largely missing because modules display context-dependent 
behavior. One cause of context-dependence is retroactivity, a phenomenon similar to loading that influences in non-trivial 
ways the dynamic performance of a module upon connection to other modules. Here, we establish an analysis framework 
for gene transcription networks that explicitly accounts for retroactivity. Specifically, a module's key properties are encoded 
by three retroactivity matrices: internal, scaling, and mixing retroactivity. All of them have a physical interpretation and can 
be computed from macroscopic parameters (dissociation constants and promoter concentrations) and from the modules' 
topology. The internal retroactivity quantifies the effect of intramodular connections on an isolated module's dynamics. The 
scaling and mixing retroactivity establish how intermodular connections change the dynamics of connected modules. Based 
on these matrices and on the dynamics of modules in isolation, we can accurately predict how loading will affect the 
behavior of an arbitrary interconnection of modules. We illustrate implications of internal, scaling, and mixing retroactivity 
on the performance of recurrent network motifs, including negative autoregulation, combinatorial regulation, two-gene 
clocks, the toggle switch, and the single-input motif. We further provide a quantitative metric that determines how robust 
the dynamic behavior of a module is to interconnection with other modules. This metric can be employed both to evaluate 
the extent of modularity of natural networks and to establish concrete design guidelines to minimize retroactivity between 
modules in synthetic systems. 
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Introduction 

The ability to accurately predict the behavior of a complex 
system from that of the composing modules has been instrumental 
to the development of engineering systems. It has been proposed 
that biological networks may have a modular organization similar 
to that of engineered systems and that core processes, or motifs, 
have been conserved through the course of evolution and across 
different contexts [1], [2], [3], [4], [5]. In addition to having 
profound consequences from an evolutionary perspective, this 
view implies that biology can be understood, just like engineering, 
in a modular fashion [6] . To predict the behavior of a network 
from that of its composing modules, it is certainly desirable that 
the salient properties of modules do not change upon connection 
with other modules. This modularity property is especially 
important in a bottom-up approach to engineer biological systems, 
in which small systems are combined to create larger ones [7] , [8] . 

Unfortunately, despite the fact that biological networks are rich 
of frequently repeated motifs, suggesting a modular organization, a 
module's behavior is often affected by its context [9]. Context- 
dependence is due to a number of different factors. These include 
unknown regulatory interactions between the module and its 
surrounding systems; various effects that the module has on the 
cell network, such as metabolic burden [10], effects on cell growth 
[11], and competition for shared resources [12]; and loading 
effects associated with known regulatory linkages between the 



module and the surrounding systems, a phenomenon known as 
retroactivity [13], [14]. As a result, our current ability of predicting 
the emergent behavior of a network from that of the composing 
modules remains limited. This inability is a central problem in 
systems biology and especially daunting for synthetic biology, in 
which circuits need to be re-designed through a lengthy and ad hoc 
process every time they are inserted in a different context [15]. 

In the phenomenon known as retroactivity, a downstream 
module perturbs the dynamic state of its upstream module in the 
process of receiving information from the latter [13], [14]. These 
effects are due to the fact that, upon interconnection, a species of 
the upstream module becomes temporarily unavailable for the 
reactions that make up the upstream module, changing the 
upstream module's dynamics. The resulting perturbations can 
have dramatic effects on the upstream module's behavior. For 
example, in experiments in gene circuits in Escherichia coli, a few 
fold ratio in gene copy number between the upstream module and 
the downstream target results in more than 40% change in the 
upstream module's response time [16]. More intriguing effects take 
place when the upstream module is a complex dynamical system 
such as an oscillator. In particular, experiments in transcriptional 
circuits in vitro showed that the frequency and amplitude of a 
clock's oscillations can be largely affected by a load [17] and 
computational studies on the genetic activator-repressor clock of 
[18] further revealed that just a few additional targets for the 
activator impose enough load to quench oscillations. Surprisingly, 
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Author Summary 

Biological modules are inherently context-dependent as 
the input/output behavior of a module often changes 
upon connection with other modules. One source of 
context-dependence is retroactivity, a loading phenome- 
non by which a downstream system affects the behavior of 
an upstream system upon interconnection. This fact 
renders it difficult to predict how modules will behave 
once connected to each other. In this paper, we propose a 
general modeling framework for gene transcription net- 
works to accurately predict how retroactivity affects the 
dynamic behavior of interconnected modules, based on 
salient physical properties of the same modules in 
isolation. We illustrate how our framework predicts 
surprising and counter-intuitive dynamic properties of 
naturally occurring network structures, which cannot be 
captured by existing models of the same dimension. We 
describe implications of our findings on the bottom-up 
approach to designing synthetic circuits, and on the top- 
down approach to identifying functional modules in 
natural networks, revealing trade-offs between robustness 
to interconnection and dynamic performance. Our frame- 
work carries substantial conceptual analogies with electri- 
cal network theory based on equivalent representations. 
We believe that the framework we have proposed, also 
based on equivalent network representations, can be 
similarly useful for the analysis and design of biological 
networks. 

adding a few targets for the repressor can restore the stable limit 
cycle [19]. Retroactivity has also been experimentally demon- 
strated in signaling networks in vitro [20] and in the MAPK cascade 
in vivo [21]. In particular, it was shown in [19] that a few fold ratio 
between the amounts of the upstream and downstream system's 
proteins can lead to more than triple the response time of the 
upstream system. 

In this paper, we provide a quantitative framework to accurately 
predict how and the extent to which retroactivity will change a 
module's temporal dynamics for general gene transcription 
networks and illustrate the implications on a number of recurrent 
network motifs. We demonstrate that the dynamic effects of 
loading due to interconnections can be fully captured by three 
retroactivity matrices. The first is the internal retroactivity, which 
accounts for loading due to intramodular connections. We 
illustrate that due to internal retroactivity, negative autoregulation 
can surprisingly slow down the temporal response of a gene as 
opposed to speeding it up, as previously reported [22]; perturba- 
tions applied at one node can lead to a response at another node 
even in the absence of a regulatory path from the first node to the 
second, having consequences relevant for network identification 
techniques (e.g., reviewed in [23]); and an oscillator design can fail 
even in the presence of small retroactivity. The other two matrices, 
which we call scaling and mixing retroactivity, account for loading 
due to intermodular connections. We illustrate that because of the 
scaling retroactivity, the switching characteristics of a genetic 
toggle switch can be substantially affected when the toggle switch is 
inserted in a multi-module system such as that proposed for 
artificial tissue homeostasis in [24]. The interplay between scaling 
and internal retroactivity plays a role in performance/robustness 
trade-offs, which we illustrate considering the single-input motif 
[5]. Using these retroactivities, we further provide a metric 
establishing the robustness of a module's behavior to interconnec- 
tion. This metric can be explicitly calculated as a function of 
measurable biochemical parameters, and it can be used both for 



evaluating the extent of modularity of natural networks and for 
designing synthetic circuits modularly. 

Our work is complementary to but different from studies 
focusing on partitioning large transcription networks into modules 
using graph-theoretic approaches [13], [25], [26]. Instead, our 
main objective is to develop a general framework to accurately 
predict both the quantitative and the qualitative behavior of 
interconnected modules from their behavior in isolation and from 
key physical properties (internal, scaling, and mixing retroactivity). 
In this sense, our approach is closer to that of disciplines in 
biochemical systems analysis, such as metabolic control analysis 
(MCA) [27], [28]. However, while MCA is primarily focused on 
steady state and near-equilibrium behavior, our approach 
considers global nonlinear dynamics evolving possibly far from 
equilibrium situations. 

This paper is organized as follows. We first introduce a general 
mechanistic model for gene transcription networks to explain the 
physical origin of retroactivity and to formulate the main question 
of the paper (System Model and Problem Formulation). We then 
provide the two main results of the paper (Results). These are 
obtained by reducing the mechanistic model through the use of 
time scale separation (leading to models of the same dimension as 
those based on Hill functions), in which only macroscopic 
parameters and protein concentrations appear. In these reduced 
models, the retroactivity matrices naturally arise, whose practical 
implications are illustrated on five different application examples. 

System Model and Problem Formulation 

We begin by introducing a standard mechanistic model for gene 
transcription networks, which includes protein production, decay, 
and reversible binding reactions between transcription factors 
(TFs) and promoter sites, required for transcriptional regulation. 
Specifically, transcription networks are usually viewed as the 
input/output interconnection of fundamental building blocks 
called transcriptional components. A transcriptional component 
takes a number of TFs as inputs, and produces a single TF as an 
output. The input TFs form complexes with promoter sites in the 
transcriptional component through reversible binding reactions to 
regulate the production of the output TF, through the process of 
gene expression (for details, see Methods). To simplify the 
notation, we treat gene expression as a one-step process, neglecting 
mRNA dynamics. This assumption is based on the fact that 
mRNA dynamics occur on a time scale much faster than protein 
production/decay [1]. In addition to this, including mRNA 
dynamics is not relevant for the study of retroactivity, and would 
yield only minor changes in our results (see Methods). 

Within a transcription network, we identify a transcriptional 
component with a node. Consequendy, a transcription network is a 
set of interconnected nodes in which node x,- represents the 
transcriptional component producing TF x,-. There is a directed 
edge from node Xj to x,- if Xj is a TF regulating the activity of the 
promoter controlling the expression of x,- [29], in which case we 
call X; a. parent of x,. Activation and repression are denoted by — > 
and H, respectively. Modules are a set of connected nodes. Modules 
communicate with each other by having TFs produced in one 
module regulate the expression of TFs produced in a different 
module. When a node x, is inside the module, we call the 
corresponding TF x,- an internal TF, while when node x,- is outside 
the module we call the corresponding TF x, an external TF. 
Further, we identify external TFs that are parents to internal TFs 
as inputs to the module. Let x, u and c denote the concentration 
vector of internal TFs, inputs and TF-promoter complexes, 
respectively. According to [30], we can write the dynamics of 
the module as 
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written as 



--N, t v(x,c,u), 



(1) 



where N sl is the stoichiometry matrix and V is the reaction flux 
vector. The reactions are either protein production/decay or 
binding/ unbinding reactions. Therefore, we partition v into r* and 
r, representing the reaction flux vectors corresponding to 
production/decay and binding/ unbinding reactions, respectively 
(see Methods). We assume that the DNA copy number is 
conserved, therefore, we can rewrite (1) as 
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where the upper left block matrix in N st is the zero matrix as DNA 
is not produced/degraded. As a result, with g(x,c) = B*r*(x,c) we 
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where the upper left block matrix is zero as DNA is assumed to be 
a conserved species. Furthermore, since r and f encapsulate the 
binding/ unbinding reactions in the module and in its context, 
respectively, the off-diagonal block matrices in the upper right 
block matrix are zero. Similarly, as r* and f* encapsulate the 
production/decay reactions in the module and its context, 
respectively, the off-diagonal block matrices in the lower left block 
matrix are zero. Finally, the stoichiometry matrix E represents 
how internal TFs of the module participate in binding/ unbinding 
reactions in the context of the module (E can be interpreted 
similarly). 

With s = Ef describing the effective rate of change of x due to 
intermodular binding reactions, we obtain 



c= Ar(x,c,u), 

x= g(x,c)+Br(x,c,u), 



(2) 



c= Ar[x,c,u), 

x= g(x,c) + Br(x,c,u) + s(x,c,u), 



(4) 



which we call the isolated dynamics of a module. 

Next, consider the case when the module is inserted into a 
network, which we call the context of the module. We represent all 
the quantities related to the context with an overbar. Let x and c 
denote the concentration vector of TFs and promoter complexes 
of the context, respectively. Furthermore, denote by r* and f the 
reaction flux vectors corresponding to production/decay and 
binding/ unbinding reactions between TFs and promoters in the 
context of the module, respectively. Then, the dynamics of the 
species in the module (c and x) and in the context (c and x) can be 



which we call the connected dynamics of a module. We refer to s as 
the retroactivity to the output of the module, encompassing retroac- 
tivity applied to the module due to the context of the module. 
Similarly, we call r the retroactivity to the input of a module, 
representing retroactivity originating inside the module. The 
general interconnection of a group of modules can be treated 
similarly (Figure 1). 

As an example of the implications of retroactivity s on the 
module's dynamic behavior, consider Figure 2. For the purpose of 
illustration, assume that i\(t) and fi(f), external inputs to Xi and 
Xi (see Methods), are periodic (in general, they can be arbitrary 




Module 1 






Module 2 


- 


r (2) 




£(3) 




Module j 



Figure 1 . The dynamics of a module depend on the module's context. Downstream modules change the dynamics of an upstream module 
by applying a load. The effect of this load is captured by the retroactivity to the output s of the upstream module, which is the weighted sum of the 
retroactivity to the input r® of the downstream modules. 
doi:10.1371/journal.pcbi.1003486.g001 
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Figure 2. The context (downstream system) affects the behavior of the module (upstream system). (A) The module in isolation. (B) The 
module in isolation displays sustained oscillations. (C) The module connected to its context. (D) Upon interconnection with its context, the dynamics 
of the module change due to the retroactivity s from its context, since some of the molecules of Xi are involved in binding reactions at node X2. As a 
result, those molecules are not available for reactions in the module, and the output of the module is severely changed. For details on the system and 
parameters, see Supporting Text SI. 
doi:1 0.1 371 /journal.pcbi.1 003486.g002 



time-varying signals). When the module is not connected to its 
context (Figure 2 A), its output is periodic (Figure 2B). Upon 
interconnection with its context (Figure 2C), due to the 
retroactivity to the output s applied by the context, the output 
of the module changes significantly (Figure 2D). Hence, connec- 
tion with the context leads to a dramatic departure of the 
dynamics of the module from its behavior in isolation. This 
example illustrates that retroactivity s significandy alters the 
dynamic behavior of modules after interconnection, therefore, it 
cannot be neglected if accurate prediction of temporal dynamics is 
required. Unfortunately, model (4) provides litde analytical insight 
into how measurable parameters and interconnection topology 
affect retroactivity. 

The aim of this paper is to provide a model that captures the 
effects of retroactivity, unlike standard regulatory network models 
of the same dimension based on Hill functions [1]. Specifically, we 
seek a model that explicitly describes the change in the dynamics 
of a module once it is arbitrarily connected to other modules in the 
network. This model is only a function of measurable biochemical 
parameters, TF concentrations, and interconnection topology. 

Results 

We first characterize the effect of intramodular connections on 
an isolated module's dynamics. We then analytically quantify the 
effects of intermodular connections on a module's behavior. 
Finally, we determine a metric of robustness to interconnection 
quantifying the extent by which the dynamics of a module are 
affected by its context. We demonstrate the use of our framework 
and its implications on network motifs taken from the literature. 

The main technical assumptions in what follows are that (a) 
there is a separation of time scale between production/ degrada- 
tion of proteins and the reversible binding reactions between TFs 
and DNA, and that (b) the corresponding quasi-steady state is 
locally exponentially stable. Assumption (a) is justified by the fact 
that gene expression is on the time scale of minutes to hours while 
binding reactions are on the second to subsecond time scale [3]. 
Assumption (b) is implicitly made any time Hill function-based 
models are used in gene regulatory networks. In addition to these 



technical assumptions, to simplify notation, we model gene 
expression as a one-step process, however, a more detailed 
description of transcription/translation would not yield any 
changes to the main results (see Methods). 

Effect of Intramodular Connections 

Here, we focus on a single module without inputs and describe 
how retroactivity among nodes, modeled by Br in (2), affects the 
module's dynamics. To this end, we provide a model that well 
approximates the isolated module dynamics, in which only 
measurable macroscopic parameters appear, such as dissociation 
constants and TF concentrations. We then present implications of 
this model for negative autoregulation, combinatorial regulation 
and the activator-repressor clock of [18]. 

Employing assumptions (a)-(b), we obtain the first main result of 
the paper as follows. Let x = (x\,X2, . . . ,x^) T denote the vector of 
concentrations of internal TFs, then the dynamics 



X=[I+R{x)\~ l h{x) 



(5) 



well approximate the dynamics of x in (2) in the isolated module 
with 



h(x)-- 



^2+H 2 {p2)-hX2 



and 
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Figure 3. Hill function and retroactivity of node \, for the most common binding types. If node \, has no parents, its node retroactivity is 
not defined. In the single parent case, node x, has one parent, y binding as an «-multimer with dissociation constant k y . In the case of independent, 
competitive and cooperative binding, node x, has two parents, y and z, binding as multimers with multimerization factors n and m, respectively, 
together with dissociation constants k v and k 2 , respectively. The total concentration of the promoter of x, is denoted by rj t . The production rates Try), 
7Cj,i, Tty and ity correspond to the promoter complexes without parents, with y only, with z only, and with both y and z, respectively. For details, see 
Supporting Text S3. 
doi:10.1371/journal.pcbi.1003486.g003 



where represents external perturbations to x, (inducer, noise, or 
disturbance, ^,-(?) = 0 unless specified otherwise), 5i is the decay 
rate of x,, and H,(p,) is the Hill function modeling the production 
rate of x,, regulated by the parents />,■ of x,. We call Ri(pi) the 
retroactivity of node x, . For the most common binding types, Figure 3 
shows the expressions of Hi(pi) and Riipi) (for their definition, see 
Methods). The binary matrix Vf has as many columns as the 
number of nodes in the module, and as many rows as the number 
of parents of x, , such that its (j,k) element is 1 if the j th parent of x,- 
is Xi, otherwise the entry is zero. That is, an entry in the following 



matrix 



xi x 2 



P«,l 
Pi,2 
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is 1 if the species indexing the corresponding row and column are 
the same, otherwise the entry is zero, yielding p, = VjX. Finally, O 
is the set of nodes having parents from inside the module. For the 
derivation of this result, see Theorem 1 in Supporting Text S2. 

We call R the internal retroactivity of the module as it describes how 
retroactivity among the nodes internal to the module affects the 
isolated module dynamics. When R = 0, we have x = h(x), the 
commonly used Hill function-based model for gene transcription 
networks [3] . It is possible to show that h(x) represents the rate of 
change of total (free and bound) TFs (see Supporting Text S2). 
Hence, (6) describes how changes in the total concentration of TFs 
h(x) relate to changes x in the concentration of free TFs. Specifically, 
to change the concentration of free TFs by one unit, the module has 
to change the total concentration of TFs by (I + R) units, as R units 
are "spent on" changing the concentration of bound TFs. Having 
R = 0 implies that the module's effort on affecting the total 
concentration of TFs is entirely spent on changing the concentration 
of free TFs. By contrast, ||i?|| — > oo implies that no matter how much 
the total concentration of TFs changes, it is not possible to achieve 
any changes in the free concentration of some of the TFs. Therefore, 
the internal retroactivity R describes how "stiff' the module is 
against changes in x due to loading applied by internal connections. 

The entries of Riipi) have the following physical interpretation. 
Consider first a module with the autoregulated node X\, that is, Xi 
has a single parent: itself. The retroactivity of node xi is Ri(x\) = a, 
where a is given in Figure 3. In this case, we obtain R(x\) = a by (6), 
so that (5) yields X\ = jxsA(jCi). Hence, the greater a, the harder to 
change the concentration of free Xi by changing its total 
concentration (the "stiffer" the node), and the temporal dynamics 
of x\ become slower. The retroactivity R\ of node Xi can be 
increased by increasing its DNA copy number r]\ or by decreasing 
the dissociation constant k\ of Xi . For a node with two parents, we 
provide the explicit formula for Rj in Figure 3 in the case of the most 
frequent binding types, so that here we simply write 



over 40% of known Escherichia coli TFs are autorepressed [29]. 
Earlier studies concluded that negative autoregulation makes the 
response of a gene faster [22]. Here, we demonstrate that in the 
case of significant retroactivity, negative autoregulation can 
actually slow down the response of a gene. To this end, consider 
a module consisting of the single node xi , and analyze first the case 
when its production is constitutive with promoter concentration 
rji, production rate constant 71^0 and decay rate <5i. Then, the 
dynamics of Xi are given by X\ —riiiii,o — SiXi. 

In the case of negative autoregulation, Xj has itself as the only 
parent. Let k\ denote the dissociation constant of Xi and assume it 
binds as a monomer repressing its own production (so that n = 1 
and 7iii=0 in Figure 3). According to Figure 3, we have 
= fl and R\ (x, ) = m together with V\ = 1 

and ® = {xi}, yielding from (6) R(x\) = R\(x\) and h(x\) = H\(xi) 
— &\X\, so that (5) results in 

* 1 = TT^(' ?i TrSk- 8i 4 (8) 

This expression indicates that negative autoregulation yields two 
changes in the dynamics. First, protein production changes from 
y\\Ti\fi to the Hill function H\(x\). Second, the dynamics are 




B Medium DNA copy number, medium production rate constant 




The diagonal entries b and e in (7) can be interpreted similarly to a, 
while the off-diagonal entries can be interpreted as follows. Having 
c> 0 means that the second parent facilitates the binding of the first, 
whereas c<0 represents blockage (d can be interpreted similarly 
with the parents having reverse roles). Therefore, we have c = d = 0 
in the case of independent binding (Figure 3), as the parents bind to 
different sites. By contrast, we have c,d < 0 in the case of competitive 
binding (Figure 3), since the parents are competing for the same 
binding sites, forcing each other to unbind. Following a similar 
reasoning, we obtain c,d>0 in the case of cooperative binding 
(Figure 3). Notice that Rj is scaled by the total concentration of 
promoter r\ h which can be changed, for example, in synthetic 
circuits by changing the plasmid copy number. 

Practical Implications of Intramodular Connections 

In order to illustrate the effects of intramodular connections, we 
consider three recurrent network motifs in gene transcription 
networks: (i) negative autoregulation of a gene, (ii) combinatorial 
regulation of a gene by two TFs, and (iii) the activator-repressor 
clock of [18]. 

Negative autoregulation. One of the most frequent network 
motifs in gene transcription networks is negative autoregulation, as 



Q High DNA copy number, low production rate constant 




0 12 3 



time [hours] 

Figure 4. Negative autoregulation can make the temporal 
response slower. Time response at a steady state fixed at x\ = 50nM. 
The red and blue plots denote the cases with and without negative 
autoregulation, respectively, whereas the green plot represents the 
case of negative autoregulation neglecting retroactivity (,ft = 0 in (8)). 
Simulation parameters are <5i = lhr~', /ci = 10nM, together with 
7ti_o = 20hr~ 1 , 7ti > o = 10hr _1 , 7ti i o = lhr _1 for A, B and C, respectively. 
To carry out a meaningful comparison between the unregulated and 
regulated systems, we compare the response time of systems with the 
same steady state. To do so, we pick the same value of \] x in the case of 
the regulated systems (i/, = 15nM, i/, =30nM, i/, =300nM for A, B and 
C, respectively), but a different one for for the unregulated system 
(1(1 =2.5nM, ;j! = 5nM, i(|=50nM for A, B and C, respectively), such 
that the steady states match (see Methods for parameter ranges). 
Decreasing m,ii (lower production rate constant) while increasing i;, 
(higher DNA copy number) results in slower response, as internal 
retroactivity increases. 
doi:10.1371/journal.pcbi.1003486.g004 
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4 6 
time [hours] 



Figure 5. Nodes can become coupled via common downstream 
targets. (A) Node x? has two parents: Xi and x%, without a regulatory 
path between them. (B) Perturbation fj applied to x\. (C) In the case of 
competitive binding, increasing the concentration of free xi yields more 
of X] bound to the promoter of X3, forcing some of the molecules of X2 
to unbind, thus increasing the free concentration x 2 . Consequently, xi 
acts as if it were an activator of X2. (D) By contrast, in the case of 
cooperative binding, when the binding of Xi must precede that of X2, 
pulses in .Ti yield pulses of the opposite sign in x 2 . Consequently, X| 
acts as if it were a repressor of x 2 . Simulation parameters are: 
1i =>/2= 10nM, >/ 3 =20nM, <5i = <5 2 = lhr -1 , 7Ci j0 =0, 7t2,o = 10hr _1 , 
k] =^2 = lnM 4 , and both Xi and x 2 bind as tetramers. 
doi:1 0.1 371 /journal.pcbi.1 003486.g005 



(1+*)" 



which is the effect of internal 



premultiplied by 
retroactivity. 

As it was shown in [22], the response time of the regulated 
system without retroactivity is smaller than that of the unregulated 
system. When considering internal retroactivity, however, the 
response time increases, as the absolute value of X\ decreases with 
increased R according to (8). Specifically, the response time with R 
is greater than without R since R>0. That is, while the Hill 
function makes the response faster, internal retroactivity has an 
antagonistic effect, so that negative autoregulation can render the 
response slower than that of the unregulated system, as illustrated 
in Figure 4. Furthermore, if 7ti,o'7i is kept constant, the response 
time of both the unregulated (blue) and the regulated system 
without retroactivity (green) remain the same, together with the 
steady states. By contrast, increasing f/j (and decreasing 7110) 
makes the internal retroactivity R greater (since R is proportional 
to JJi), while the contribution of the Hill function remains 
unchanged. As a result, the response of the regulated system with 
retroactivity (red) becomes slower as we increase rj l (and decrease 
7ti,o). This is illustrated in Figure 4 with different (rji,ni^) pairs. 
Note that 7ti o can be decreased, for example, by decreasing the 
ribosome binding site (RBS) strength, whereas r]\ can be increased 
by increasing the gene copy number. 

Combinatorial regulation. As a second example, we 
consider a single gene co-regulated by two TFs (Figure 5A). This 



topology appears in recurrent network motifs, such as the 
feedforward-loop, the bi-fan and the dense overlapping regulon 
[5] . Here, we show that a perturbation introduced in one of the 
parents (blue in Figure 5A) can affect the concentration of the 
other parent (red node in Figure 5A), even in the absence of a 
regulatory path between the two. 

Referring to (5)-(6), note that X3 is the only node with parents 
(3> = {x 3 }), so that R(x) = V[R 3 V 3 . Using (6) with 



1 0 0 
0 1 0 



and R 3 (xi,Xz)-- 



b c 
d e 



where the entries of Rt, are given in Figure 3 (depending on the 
binding type at X3) together with H-}{x\,xi), the dynamics in (5) 
take the form 



x 2 

w 



1+e 
(l+b)(l+e)-cd 

d 

(\+b)(\+e)-cd 

0 



0 



(l+6)(l+e)-crf 
Lb* 0 

0 1 



-1 



[I + R(x)\ 

I Ci+7ii,o>/i-6ijci \ 

ft2,0'?2- S 2*2 
V ff 3 (*l,-Y2)-8 3 *3 ) 

HA 



This expression implies that unless d = 0, a perturbation 
(Figure 5B) in X\ yields a subsequent perturbation in x%. In the 
case of independent binding, we have d = 0, and as a result, no 
perturbation is observed in Xj. In the case of competitive binding, 
instead, we have d<0, so that perturbations £j in X\ yield 
perturbations of the same sign in Xi, that is, Xi acts as if it were an 
activator of X2 (Figure 5C). In the case of cooperative binding, 
instead, we have d>0. As a result, perturbations in X\ yield 
perturbations in X2 of opposite sign (Figure 5D), which implies that 
Xi behaves as if it were a repressor of X2. As d is proportional to 173 
(Figure 3), higher DNA copy number for X3 yields greater pulses in 
X2 subsequent to an equal perturbation in X\. Interestingly, if we 
view X2 as the output of the module, the module has the 
adaptation property with respect to its input X\ (or i^). That is, 
retroactivity enables to respond to sudden changes in input stimuli, 
while adapting to constant stimulus values. 

Activator-repressor clock. One common clock design is 
based on two TFs, one of which is an activator and the other is a 
repressor [18], [31], [32]. Here, we illustrate the effect of internal 
retroactivity on the functioning of the clock design of [18] depicted 
in Figure 6A. In particular, X] activates the production of both 
TFs, whereas X2 represses the production of xi through 
competitive binding. Consequently, the network topology is 
captured by the binary matrices V\=l and Vi = [ 1 0 ] , whereas 
h(x) and R(x) can be constructed by considering (Hi(xi,X2), 
Hi(x\)) and (R\{x\,X2), Rjixi)), respectively, in Figure 3. Here, 
we write R2 = a, while the entries of R\ are denoted by b, c, d and 
e, as in (7). Then, we obtain that (5) takes the form 
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Figure 6. Neglecting internal retroactivity could falsely predict 
that the activator-repressor clock will display sustained 
oscillations. (A) The module consists of the activator protein xi 
(dimer) and the repressor protein x 2 (tetramer), with dissociation 
constants k\ and k%, respectively. (B) An extra node X3 is introduced as 
a target for the repressor. (C) Without accounting for internal 
retroactivity, the module in A exhibits sustained oscillations. (D) When 
internal retroactivity is included for the module in A, however, the 
equilibrium point is stabilized and the limit cycle disappears. (E) 
Oscillations can be restored by applying a load on the repressor (module 
in B) with concentration ;/, so that the repressor dynamics are slowed 
down. Simulation parameters: jti.o = 0.04hr _1 , jri.i =2500hr~\ 
7ti,2=0hr- 1 , 712.0 = 0.004^-', n 2 .i = lOChr" 1 , <5, = lhr" 1 , fe = 0.6hr- 1 , 
1i =Jj2=2nM,fci =lnM 2 , k 2 = InM 4 and ;; = 40nM. 
doi:1 0.1 371 /journal.pcbi.1 003486.g006 



l+e c 

(l+a + b)(l+e)-cd (l+a+A)(l+e)-cd 

d l+a+b 

(l+a+b)(l+e)-cd (1 + a+b)(l +e)-cd 



[I + R(x)\ 



-1 



(9) 



Hi(x u x 2 )-5ixi 
H 2 (x l )-8 2 X2 



It was previously shown [33] that the principle for the clock to 
oscillate is that the activator dynamics are sufficiently faster than 
the repressor dynamics (so that the unique equilibrium point is 
unstable). Equation (9) shows that the activator and repressor 
dynamics are slowed down asymmetrically (diagonal terms in 

[I + R(x)] 1 ), and that they become coupled (off-diagonal terms in 
[7+J?(^c)] , c,d^=Q), because of internal retroactivity. In 
particular, in the case when c,d« 1 +e« 1 +a + b, the activator 
would slow down compared to the repressor. Based on the 
principle of functioning of the clock, we should expect that this 
could stabilize the equilibrium point, quenching the oscillations as 
a consequence. In fact, oscillations disappear even if the circuit is 
assembled on DNA with a single copy (f/j = t] 2 = 2nM), as it can be 



observed in Figure 6D. Therefore, accounting for internal 
retroactivity is particularly important in synthetic biology during 
the design process when circuit parameters and parts are chosen 
for obtaining the desired behavior. An effective way to restore the 
limit cycle in the clock, yielding sustained oscillations, is to render 
the repressor dynamics slower with respect to the activator 
dynamics. This can be obtained by adding extra DNA binding 
sites for the repressor [19], as shown in Figure 6B. In fact, in this 
case, we have R^(x 2 )>0 given in Figure 3, which, due to (5), will 
yield the following change in (9): instead of e, we will have 
e + R3 > e, rendering the dynamics of the repressor slower with 
respect to the activator dynamics. As a result, the equilibrium 
point becomes unstable, restoring the limit cycle, verified by 
simulation in Figure 6E. Further studies on specific systems have 
investigated the effects of TF/ promoter binding on the dynamics 
of loop oscillators, such as the repressilator [34] . 

Effect of Intermodular Connections 

After investigating how retroactivity due to intramodular 
connections affect a single module's dynamics, we next determine 
how the dynamics of a module change when the module is inserted 
into its context. To this end, we first extend the model in (5) to the 
case in which the module has external TFs as inputs. Hence, let 
u = (u\,u 2 , . . . ,u w) T denote the concentration vector of TFs 
external to the module. With this, we obtain that the dynamics 

x=[I + R(x,u)] ~~ 1 [h{x,u) - Q(x,u)u] =f(x,u,u) (10) 
well approximate the dynamics of x in (2) with 



f E VfR i ( Pi )D i 
Q(x,u)=< {i|x,- 6 (*nn)} 

I OjVx w 



if n#0, 



if n= 



an 



where Q is the set of nodes having parents from outside the 
module (external TFs), and the binary matrix D, has as many 
columns as the number of inputs of the module, and as many rows 
as the number of parents of x, , such that its (j,k) element is 1 if the 
j th parent of x,- is u^, otherwise the entry is zero. That is, an entry 
in the following matrix 



Ui 112 



D, 



P/,1 

P/.2 



is 1 if the species indexing the corresponding row and column are 
the same, otherwise the entry is zero, yielding p, = [ Vj Z), ] 
(x T u J ) . Furthermore, note that in the presence of input u, 
both h(.) and R(.) given in (6) depend on x and u, as some of the 
parents of internal TFs are external TFs. For the derivation of this 
result, see Theorem 2 in Supporting Text S2. 

Before stating the main result of this section, we first provide the 
interpretation of Q. Recall that h(x,u) = 0 implies that the total 
concentrations of internal TFs are constant. In this case, (10) 
reduces to x = — (I + R) ~ 1 Qii, where x is the concentration 
vector of free internal TFs. This means that the concentrations of 
free internal TFs can still be changed subsequent to changes in the 
external TFs (input), despite the fact that the total concentration 
(free and bound) of internal TFs remains unaffected. Therefore, Q 
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captures the phenomenon by which external TFs force internal 
TFs to bind/unbind, for instance, by competing for the same 
binding sites. Having Q = 0 means that external TFs do not affect 
the binding/unbinding of internal TFs, which is the case, for 
example, when all bindings are independent. Thus, we call Q the 
external retroactivity of the module. 

The main result of this section describes how the context of a 
module affects the module's dynamics due to retroactivity. 
Specifically, we consider the module of interest and we represent 
the rest of the network, the module's context, as a different 
module. As previously, we use the overbar to denote that a 
quantity belongs to the context. With this, we obtain that the 
dynamics 



interconnection become 



I+(I+R)- l S 

{i+Ry l M 

f(x,Ux,Ux) 
f(x,Ux,Ux) 



(i+Ry l M 

I + il + R)-^ 



(12) 



isolated dynamics 
of the module and 
of its context 

well approximate the dynamics of x and x in (4) in the module 
connected to the context with 



J" £ [DtU? Rfp t )DtU if a *0, 

S(x,x)= < {i|x,-eQ} 

I <W if 0 = 0, 

£ [DiU] T RtiptWu if*rW0, 



M(X,X)= < {i|x,e(Opl£2)} 



if onn=0, 



where N and N denote the number of nodes in the module and in 
the context, respectively. Furthermore, the binary matrix U has as 
many rows as the number of inputs of the module, and as many 
columns as the number of nodes in the context, such that its (J,k) 
element is 1 if the j th input of the module is the k th internal TF of 
the context (u, = X/t), otherwise the entry is zero. That is, an entry 
in the following matrix 



Xi x 2 



U-- 



ui 

"2 



is 1 if the species indexing the corresponding row and column are 
the same, otherwise the entry is zero, yielding u=Ux. The 
quantities corresponding to the context, that is, S, M and U are 
defined similarly with the only difference that variables with and 
without overbar have to be swapped (for instance, N and N have 
to be swapped in (13)). For the derivation of this result, see 
Theorem 3 in Supporting Text S2. 

We next provide the interpretation of the scaling and mixing 
retroactivity. The reduced order model (12) describes how 
retroactivity between the module and the context affects each 
other's dynamics. Note that zero matrices S, M, S and M lead to 
no alteration in the dynamics upon interconnection. To further 
deepen the implications of these matrices and their physical 
meaning, note that when M = 0, the dynamics of the module after 



x= [l+{I + R)- l S~\ l f(x,Ux,Ux), 



(14) 



isolated dynamics 
of the module 



that is, S determines how the isolated dynamics of the module get 
"scaled" upon interconnection. Therefore, we call 5* the scaling 
retroactivity of the context, accounting for the loading that the 
context applies on the module as some of the TFs of the module 
are taken up by promoter complexes in its context (we obtain S = 0 
if nodes in the context do not have parents in the module, that is, if 
£2 = 0). Since the dynamics of the context enter into the module's 
dynamics through M, we call M the mixing retroactivity of the 
context, referring to the "mixing" of the dynamics of the module 
and that of its context. The mixing retroactivity M establishes how 
internal TFs force external TFs to bind/ unbind, so that M = 0 can 
be ensured if the binding of parents from the module is 
independent from that of the parents from the context. This 
holds if nodes in the context are not allowed to have parents in 
both the module and in the context (£2n© = 0). When M #0, a 
perturbation applied in the context can result in a response in the 
upstream module, even without TFs in the context regulating TFs 
in the module, leading to a counter-intuitive transmission of signals 
from downstream (context) to upstream (module). 

With this, we can explain the simulation results in Figure 2D by 
analyzing S and M for the system in Figure 2C. Let R\=a and let 
Rj be defined as in (7), where a, b, c, d and e are given in Figure 3. 
Then, we have R = a by (13) and S = b and M = (c 0) by (6). 
Hence, expression (12) yields 



1 



x\ ■■ 



-/(*)- 



z fi(x,x,x), 



(15) 



describing the dynamics of X\ upon interconnection with its 
context, where x\ =f(x) and x\ =f\{x,x,x) describe the dynamics 
of x\ and x,\, respectively, when the module and the context are 
not connected to each other. If M = 0, then (15) reduces to 



x\ ■■ 



l+a 



1 



/(*). 



that is, the context rescales the dynamics of the module. The 
smaller (1 +a)/(l +a + b), that is, the greater the scaling 
retroactivity b = S, the greater the effect of this scaling. Note that 
since the scaling factor is smaller than 1 (unless the scaling 
retroactivity is zero, i.e., b = 0), the effect of the scaling 
retroactivity of the context in this case is to make the temporal 
dynamics of the module slower. 

Once Mj^O, in addition to this sclaing effect, the dynamics of 
the context appear in the dynamics of the module (Figure 2D). 
Referring to (15), we can quantify the effect of the context on the 
module, considering the ratio c/(l+a). The greater c/(l+a), 
that is, the greater M, the stronger the contribution of the context 
compared to that of the module to the dynamics of the module 
upon interconnection. Here, both b and c increase, for instance, 
with the copy number of X2 (Figure 3). 

Connecting the module to its context such that Xj and Xi are 
competing for the same binding sites is less desirable than 
employing independent binding, as the dynamics of the context 
(downstream system) can suppress the dynamics of the module 
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(upstream system). Dismantling the dynamics of the module will 
"misinform" other downstream systems in the network that are 
regulated by Xi . From a design perspective, multi-module systems 
should be designed and analyzed such that the modules have zero 
mixing retroactivity. This can be achieved, for instance, by 
avoiding non-independent binding at the interface nodes (at x 2 m 
Figure 2B). However, since completely independent binding can 
be hard to realize in the case of combinatorial regulation, nodes 
integrating signals from different modules should not be placed 
into the input layer (nodes having parents from other modules), 
yielding QC\Q> = 0. This can be achieved by introducing an extra 
input node in the downstream module (see Supporting Figure SI). 

Next, we quantify the difference between the isolated and 
connected module behavior. In particular, we provide a metric of 
the change in the dynamics of a module upon interconnection 
with its context, dependent on R and S and under the assumption 
that M = 0 (obtained, for instance, by avoiding mixing parents 
from the module and the context). The isolated dynamics of the 
module can be well approximated by the reduced order model 
x=f{x,u,u) in (10), and let x{i) denote its solution. Once we 
connect the module to its context, its dynamics change according 
to (14), which we write as x=f(x,u,u) and let x(t) denote the 
corresponding solution. Using the sub-multiplicative property of 
the induced 2-norm, we have that the percentage change of the 
dynamics upon interconnection can be bounded from above as 
follows: 



l ^'"'" ) - /Mll2 <^)=||[/ + (/ + ^)- 1 S]- I -/||,(16) 



\[f(x,u,ii)\\ 2 



Furthermore, with ji>0 independent of x and u, such that 
fi(x,u) < fi, we obtain that 



\\x(t)-x(t)\\ 2 = 0(fi), 



partitionings based on other measures (e.g., edge betweenness 
[35], its extension to directed graphs with nonuniform weights 
[36], round trip distance [25] or retroactivity [37]) with respect to 
robustness to interconnection. From a design point of view, one 
can design multi-module systems such that internal, scaling and 
mixing retroactivities allow for low values of fi, leading to modules 
that behave almost the same when connected or isolated. 

Practical Implications of Intermodular Connections 

We next illustrate the effect of intermodular connections on the 
dynamics of interconnected modules, considering both a synthetic 
genetic module that is being employed in a number of applications 
and a natural recurring network motif. 

Toggle switch. Here, we consider the toggle switch of [38], a 
bistable system that can be permanently switched between two 
steady states upon presentation of a transient input perturbation. 
This module has been proposed for synthetic biology applications 
in biosensing (see, for example, [34], [39]). In this paper, we 
consider the toggle switch inserted into the context of the synthetic 
circuit for controlling tissue homeostasis as proposed in [24], and 
investigate how the context of the toggle affects its switching 
characteristics. Figure 7A illustrates the toggle switch in isolation, 
whereas Figure 7B shows the configuration when connected to the 
context [24]. Note that all nodes, both in the toggle switch and in 
its context, have a single parent. Therefore, H\{x2), H 2 (x\), 
H\(x\), H 2 (x\), H 3 (x 2 ), and similarly, R\(x 2 ), R 2 (xi), Ri(xi), 
R 2 (x\), R-3(x 2 ) are given in Figure 3. 

We first consider the model of the toggle switch when not 
connected to its context (Figure 7A). Since the toggle switch has no 
input, its isolated dynamics are described by (5), where 
0 = {xi,x 2 } ; Fi = [0 l]andF 2 = [l 0]yield 



l 

0 



\+R x (x 2 ) 



[I+R(x)Y 



Ci +H l (x 2 )-d i xi 
H 2 (xi)-8 2 x 2 



=/(*)■ 



that is, fi provides an upper bound on the percentage change in 
the dynamics of the module, and on the difference in the 
trajectories of the module upon interconnection with its context. 
Furthermore, by using the properties of the induced 2-norm, we 
obtain that we can pick 



fi = max =- 

x ' x ffminU +R) — ffmax(S) 



provided that a mdx (S)<a m i n (I + R) for all x and x, where 
Cmin(7 + -R) an d ffmax(S) denote the smallest singular value of 
(I + R) and largest singular value of S, respectively. For the 
mathematical derivations, see Supporting Text S2. This suggests 
that the module becomes more robust to interconnection by 
increasing rnin Y - a m { n (I + R) or by decreasing max x - a mdx (S). 

Such a metric can be used both in the analysis and in the design 
of complex gene transcription networks as follows. Given any 
network and a desired module size N (number of nodes within the 
module), we can identify the module that has the least value of p., 
that is, the module with the greatest guaranteed robustness to 
interconnection. Furthermore, we can also evaluate existing 



Next, we consider the toggle switch connected to its context 
(Figure 7B). As nodes in the toggle switch have no parents from 
outside it, we have S = 0 and M = 0 by (13). Nodes in the context 
have no parents in the context, leading to _R = 0 from (6), and to 
2 = 0, referring to (11). With this, the isolated dynamics of the 
context are given by 









-I 







Hi(xi)-5\xi 
H 2 (xi)-5 2 x 2 

H3(X 2 )-$3X3 , 
h(x,u) 



=f(x,u), 



according to (10). The fact that nodes in the context do not have 
mixed parents from the toggle switch and from the context results 
in M = 0 from (13). With Q = {xi,x 2 ,x 3 }, D X =D 2 = [\ 0], 
Z> 3 = [0 1], and U = I we obtain 



Ri+R 2 
0 



0 

R 3 
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time [hours] 

Figure 7. Effects of the context on the switching characteristics of the toggle switch. (A) The toggle switch in isolation. (B) The toggle 
switch connected to its context [40], (C) A narrow pulse in £j (input perturbation in x\, depicted in green) causes the isolated toggle to switch 
between the two stable equilibria. (D) When connected to the context, the same pulse is insufficient to yield a switch. (E) With a wider pulse, the 
switching is restored (however, dynamics are slower compared to the isolated case). Simulation parameters: both Xi and x? bind as dimers, 
>h = '?2 = 10nM, k\ =1(2 = InM 2 , (5i =<?2 = lhr -1 , 711,0 = 712,0 = 10hr _1 , fjj =fj 2 = )j 3 = 5nM, and the height of the input perturbation pulse is 
^ = lOnMhr" 1 . 

doi:1 0.1 371 /journal.pcbi.1 003486.g007 



As a result, with a(x i ) = (R l + R 2 )/(l + R 2 )>0 and [l(x 2 ) = R 3 / 
(1 + > 0, the dynamics of the toggle switch once connected to 
the context (Figure 7B) are given by 



" 1 


0 ~ 




( 


0 


1 




1+/L 





Si +H l (x 2 )-S l x 1 
l + R 2 
H 2 (x\)-& 2 x 2 
l+Rl 



v 71 , s v ' 

+ + f(x) 

according to (12), so that the dynamics of X\ and x 2 are unaffected 
if a = 0 and /? = 0, respectively. When a,/?>0, both the x\ and x 2 
dynamics become slower upon interconnection, so that the 
response to an input stimulation will also be slower. As a 
consequence, upon removal of the stimulation, the displacement in 
the toggle state may not be sufficient to trigger a switch. This is 
illustrated in Figure 7C-D. In order to recover the switch, a wider 
pulse is required (Figure 7E) to compensate for the slow-down due 
to the context (also, note that the switching dynamics are slower 
than in the isolated case). As a result, even if the toggle had been 
characterized in isolation, it would fail to function as expected 
when inserted into its context. Note that we have ji = max^^ 
a(xi) P(x 2 ) 

(- ; — -,- — — -), where a represents the amount of load on 

l + a(xi)'l + £(x 2 ) 
Xi imposed by the context compared to that by the module, and ft 
can be interpreted similarly. The greater a (or /J), the slower the 
dynamics of X\ (of xi) become upon interconnection with the 
context. Greater a and ji yield greater p., suggesting decreased 
robustness to interconnection, verified by the simulation results. 

Single-input motif. As a second example, we focus on a 
recurrent motif in gene transcription networks, called the single- 
input motif [5]. The single-input motif is defined by a set of 
operons (context) controlled by a single TF (module), which is 
usually autoregulated (Figure 8A). It is found in a number of 
instances, including the temporal program controlling protein 
assembly in the flagella biosynthesis [40] . Here, we show that the 
dynamic performance (speed) of the module and its robustness to 



interconnection with its context are not independent, and that this 
trade-off can be analyzed by focusing on the interplay between the 
internal retroactivity R of the module and the scaling retroactivity 
S of the context. 

The isolated dynamics of the module are given in (8), which we 
write here as Xi=f(x\). Furthermore, we have Di = \ for 
!=1,2,...,/ and (7=1, so that S(x 1 )= ]T)' =1 Rfat) by (13), 
where / is the number of nodes in the context and Ri(x\) is given 
in Figure 3 (single parent). Consequently, upon interconnection, 
the dynamics of the module change according to (14) as 

l+i?(xi ) + S(xi) v 

N v effect of 

effect of the context the context 

where /J,(x{) = S(x{)/[l + R(xi) + S(xi)] and equals the expres- 
sion in (16). The smaller fl(xi), the more robust the module to 
interconnection. Note that R(xi) is proportional to rj\, therefore, 
while increasing R(x\) makes the module slower (Figure 8B), it 
also makes it more robust to interconnection (Figure 8C). It was 
previously shown that negative autoregulation increases robustness 
to perturbations [41]. Here, we have further shown that increasing 
the internal retroactivity R of the module provides an additional 
mechanism to increase robustness to interconnection, at the price 
of slower response. For a fixed steady state (the product of n\fi and 
J7 j is held constant), smaller tt^o yields greater r\ x , that is, increased 
R and, in turn, smaller fi. From a design perspective, if speed is a 
priority, one should choose a strong RBS with a low copy number 
plasmid, or alternatively, a promoter with high dissociation 
constant k\. By contrast, if robustness to interconnection is 
central, a weak RBS with a high copy number plasmid (or with 
low k\) is a better choice. If both speed and robustness to 
interconnection are desired, other design approaches may be 
required, such as the incorporation of insulator devices, as 
proposed in other works [42]. 
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Figure 8. Internal retroactivity makes a module more robust to 
interconnection at the expense of speed. (A) The module consists 
of a single negatively autoregulated node, whereas the context 
comprises / nodes repressed by the TF in the module. (B) The internal 
retroactivity R of the module increases with the DNA copy number v\ x 
of xi. As a result, the module becomes slower as R increases. (C) The 
percentage increase in the response time of the module decreases with 
r\ x , that is, internal retroactivity R increases the robustness to 
interconnection. Simulation parameters: 5i = lhr _1 , /:i = 10nM and 
tti_o is changed such that x\ =50nM at the steady state. The context 
contains / nodes each with DNA concentration r\j = InM for i= 1,2, . . . ,/ 
(low load: /=10; medium load: 1 = 20; high load: / = 50). The response 
time is calculated as the time required to reach 50% of the steady state 
value. 

doi:1 0.1 371 /journal.pcbi.1 003486.g008 



Remark. The above presented trade-off between robustness 
to interconnection and dynamic performance can be observed also 
in electrical systems. To illustrate this, consider the electrical 
circuit in Figure 9A consisting of the series interconnection of a 
voltage source /, a resistor R and a capacitor C, in which the 
output voltage is W. The speed of the circuit can be characterized 
by its time constant t = RC: the greater T, the slower the response. 
Upon interconnection with its context, represented by the 
capacitor C, the time constant of the system changes to 
T = R(C+C), while the steady state remains the same. Note that 
the percentage change in x decreases with C, making the module 
more robust to interconnection, at the expense of slower response 
when isolated. 

To further generalize the analogy between electrical systems 
and gene transcription networks [43], consider the electrical 
circuits in Figure 9B. When the module is not connected to its 
context, we have w=f and W=f, which changes to 





Module 



Context 



B 




Module 



Context 



Figure 9. Analogy with electrical systems. (A) The module consists 
of the series interconnection of a voltage source /, a resistor R and a 
capacitor C. The speed of the module can be characterized by the time 
constant %=RC, which increases upon interconnection with the 
context. The greater C, the slower the module in isolation, but the 
smaller the percentage change in its speed upon interconnection. (B) 
According to the fundamental theorem by Thevenin [48], any linear 
electrical network can be equivalently represented by a series 
interconnection of a voltage source and an impedance. As a result, a 
generic module consists of the series interconnection of a voltage 
source / and an impedance Z, and similarly, any context can be 
represented with the series interconnection of / and Z. 
doi:10.1371/journal.pcbi.1003486.g009 



1 1 

1+Z/Z 1+Z/Z 

1 1 



1+Z/Z 1+Z/Z, 



upon interconnection. This relationship is conceptually analogous 
to (12). That is, the module is robust to interconnection with its 
context if Z is small compared to Z, whereas the genetic module is 
robust to interconnection with its context if S and M are "small" 
compared to R. Therefore, R is conceptually analogous to 1/Z 
(output admittance), whereas S and M play a role similar to 1 /Z 
(input admittance). 

Discussion 

In this paper, we have focused on retroactivity, one source of 
context-dependence, and demonstrated that the internal, scaling, 
and mixing retroactivity provide missing knowledge that captures 
loading effects due to intramodular and intermodular connections. 
The internal retroactivity quantifies the effect of intramodular load, 
applied by nodes within a module onto each other because of 
binding to promoter sites within the module. Given a module of 
interest, the effects of intramodular loading on the module's 
dynamics are captured by equations (5)-(6), in which one needs to 
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replace the specific expressions of the Hill functions and node 

retroactivities Ri(p,) provided in Figure 3, and the binary matrices 
Vj encoding the network topology. The scaling retroactivity 
accounts for the intermodular loading that the context applies on 
a module, due to having some TFs of the module bound to 
promoter sites belonging to the context. The mixing retroactivity 
couples the dynamics of the module and that of the context upon 
interconnection, and it is non-zero when TFs from different 
modules bind non-independently at promoter sites. The effects of 
intermodular loading are captured by equations (10) — (1 3). To 
obtain this description, it is sufficient to consider the Hill functions 
Hi(pi) and node retroactivities RApi) provided in Figure 3, together 
with the binary matrices K,, Z>, and U representing the network 
topology. In general, the effects of the retroactivity matrices tend to 
increase with increased DNA copy number and/ or decreased 
dissociation constants. 

We have illustrated that accounting for retroactivity reveals 
surprising dynamical properties of modules and, at the same time, 
can aid design. For example, negative autoregulation, depending 
on the gene copy number and production rates, can slow down the 
response of a system instead of speeding it up. A gene can respond 
to a perturbation applied to a different gene even in the absence of 
a regulatory path between the two genes. We have shown that this 
can occur when a group of TFs co-regulate common targets and 
these common targets are found in abundance. This type of motif, 
referred to as the dense overlapping regulon, is highly frequent in 
natural regulatory networks [1]. As a result, system identification 
techniques based on perturbation analysis [23] could erroneously 
identify non-existent regulatory linkages if retroactivity is not 
accounted for in the corresponding models. An activator-repressor 
clock on low copy DNA plasmid displays sustained oscillations 
when internal retroactivity is neglected, while oscillations are 
quenched once internal retroactivity is accounted for. However, by 
carefully adjusting the module's internal retroactivity through the 
addition of DNA load for the repressor, we can restore oscillations. 
A genetic toggle switch that can be flipped by a transient external 
stimulation requires a substantially longer stimulation to be flipped 
once it is connected to just few downstream targets. These facts are 
relevant, in particular, when designing synthetic biology circuits 
and multi-module systems. 

Similar to synthetic systems, natural systems are subject to 
retroactivity. For example, clocks responsible for circadian 
rhythms have a large number of downstream targets [44], [45], 
which, in turn, may apply substantial load. This load can affect the 
amplitude and frequency of oscillations of the clock as well as the 
stability of the corresponding limit cycle. This suggests that natural 
systems may have evolved to use retroactivity in advantageous 
ways such as using it to properly tune the dynamic behavior of a 
module without changing the module's components. This 
hypothesis is further supported by the fact that there are a large 
number of TF binding sites on the chromosome that do not have a 
regulatory function [46], [47]. These sites have an impact on the 
temporal response of TFs, and could therefore be exploited by 
nature to further control the dynamics of gene regulation. More 
generally, retroactivity provides means for information to travel 
from downstream targets to upstream regulators, therefore 
establishing indirect connections. In highly interconnected topol- 
ogies, this information transfer can result in previously unknown 
ways of realizing sophisticated functions. One such example that 
we have provided is the adaptation function that topologies such as 
the dense overlapping regulon can realize by virtue of having 
nodes co-regulate multiple downstream targets. 

Based on the three retroactivity matrices, we provided a metric 
of robustness to interconnection, quantifying the percent change 



between the dynamics of a module in isolation and once connected 
to other modules. This metric is an explicit function of measurable 
parameters and becomes smaller when a module's internal 
retroactivity is large compared to the scaling retroactivity of the 
modules it connects to. This interplay may help uncover trade-offs 
in natural systems, providing a new angle for understanding 
natural principles of network organization. From an engineering 
perspective, we have provided quantitative design tools that can be 
employed in synthetic biology to appropriately match the internal 
and scaling retroactivity of connected circuits to preserve the 
circuits' behavior upon interconnection with different contexts. 
Our metric of robustness to interconnection further allows to 
evaluate the extent of modularity of a dynamical module, possibly 
enabling the discovery of previously unknown core processes. Our 
metric could be employed by currently available methods for 
partitioning networks into modules. Specifically, to evaluate 
connectivity, these methods rely on several metrics, for instance, 
edge betweenness [35], its extension to directed graphs with 
nonuniform weights [36], round trip distance [25] or retroactivity 
[37]. The metric of robustness to interconnection that we have 
introduced can enhance these methods by providing a way to 
evaluate modules on the basis of their functional robustness to 
interconnection in addition to distinguishing them at the 
connectivity level. 

The framework that we have proposed carries substantial 
conceptual analogies with the electrical circuit theory established 
by Thevenin [48], which has been used for more than a hundred 
years to analyze and to design electrical networks. Within this 
theory, each circuit has an equivalent input and output impedance 
(conceptually analogous to the scaling/mixing and internal 
retroactivity, respectively), and an equivalent energy source 
(playing a role similar to the isolated module dynamics). This 
theory has been instrumental for answering key questions in the 
analysis and design of electrical networks including, for example, 
how the output of a circuit changes after it is interconnected in a 
network; how to design circuits to maximize the power transfer 
upon connection (impedance matching); and how to design 
circuits whose input/ output response is unaffected by loads. We 
believe that the framework proposed in this paper can be used in a 
similar way for the analysis and design of gene regulatory 
networks. 

Although our framework can be applied to a general gene 
transcription network, there are a number of aspects that it does 
not currendy capture. These include post-translational protein 
modifications, such as phosphorylation, which are present in many 
regulatory networks and may potentially affect retroactivity. 
Including these will require to extend our framework to mixed 
gene transcription and signaling network models, leading to 
systems with multiple time scales. Furthermore, the transcription 
and translation processes use shared resources such as RNA 
polymerase and ribosomes, which may create couplings among 
unconnected circuits [17]. The dynamics of shared resources has 
not been included in our modeling framework and will be the 
focus of our future work. 

Methods 

Detailed Description of the System Model 

The production of TF x, is regulated by its parents p, i>P!'2! • ■ - : 
they bind to the promoter of x,-, and form complexes C,,i, Cy, . . . 
with the promoter. Each of these complexes, in turn, produce x,- 
with a different rate, where we use a one-step production process 
encapsulating both transcription and translation [3]. As a result, 
the reactions we consider for node x, are 
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Cy-^cy+x,-, (17) 



modeling the following physical phenomena. We denote by 5, 
protein decay, whereas represents the production rate that 
may be due to external inputs or perturbations (inducer, noise or 
disturbance). The second reversible reaction in (17) describes the 
binding of parent p, / with multimerization factor Wy to promoter 
complex Cy forming complex Cy-, where ay^ and Py^j are the 
association and dissociation rate constants, respectively. Further- 
more, each promoter complex c,y will contribute to the production 
of x, through the production rate constant Tty, modeled by the 
third reaction in (17). This production rate constant is a lumped 
parameter that incorporates features such as the RBS strength and 
the promoter strength. Finally, we assume that the total 
concentration of the promoter, denoted by rj h for each transcrip- 

tion component is conserved, so that 17, = z2jLo c 'J> w l lere C, is 
the number of possible complexes formed with the promoter of X;. 
This concentration is proportional to the concentration of copies 
of the promoter, which can be controlled, for example, by 
changing plasmid copy numbers in synthetic systems. 

The reaction flux vector v contains all the reactions in the 
system, that is, binding/ unbinding and protein production/ decay. 
Given that binding/unbinding reactions occur on a much faster 
time-scale than protein production/decay [1], we partition v into 
r* and r, where r* contains the slow processes, whereas r is 
composed of the fast reactions, that is, 



\ 
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(18) 



approximately 1 — 5nM concentration. By [52], a typical value of 
the dissociation constant of bacterial promoters is k=lnM, 
whereas [22] suggests /:=10nM, and experimentally obtained 
values are provided in [53]. One of the most widely used high copy 
number vectors is the pUC plasmid [54], which can have 
hundreds of copies per cell [55] . A frequendy used medium copy 
number plasmid is pl5A with a few dozen copies per cell [56], 
whereas pSC 1 0 1 is regarded as a low copy number plasmid with 
only a few copies per cell [56]. Finally, since the lifetime of a 
protein is on the order of a cell-cycle [1], we have 
<5 = 0 . 3 — 1 . 2hr 1 [50]. The typical range of macroscopic param- 
eters in Escherichia coli is summarized in Table 1 . 

If we had not neglected mRNA dynamics, there would be three 
different time scales in the system. Binding and unbinding 
reactions occur on the time scale of seconds (or even subseconds) 
[1], representing the fastest time scale. The intermediate time scale 
is that of mRNA dynamics, as the average lifetime of mRNA is on 
the time scale of minutes [57], [58], [59]. Finally, the dynamics of 
proteins evolve on the slowest time scale (hours). As we are 
interested in describing the dynamics of the system on the time 
scale of gene expression, the concentration of promoter complexes 
and mRNA transcripts can be both approximated with their quasi- 
steady state values, leading to the models we have proposed in the 
paper. However, we would like to point out that including mRNA 
dynamics would not change anything substantial in the results and 
it would simply add N more ODEs to the ODE model of an 
A-node module without any effects on the retroactivity matrices 
(shown in [8] considering a specific example). 

Definition of H t {pi) and R/ip/) 

First, note that A in (2) has a block diagonal structure, yielding 



\c N J 



A x 0 
0 A 2 

0 0 



ri{p2,ci) 



(19) 
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Biochemical Parameters 

Since the production of a typical protein takes approximately 
5 minutes [1], and a few dozen mRNAs can be transcribed from 
the same gene simultaneously by [49] , and similarly, a few dozen 
proteins can be translated from the same mRNA at the same time 
by [50], the effective production rate of protein from a gene can be 
as high as n k lOOOOhr -1 . This value can be arbitrarily decreased, 
for instance, by decreasing the RBS strength in synthetic circuits. 
The cell volume of Escherichia coli is typically between 
0.34— 1.32um 3 by [51], so that 1 molecule/cell corresponds to 



where r,(pi,Cj) denotes the reaction flux vector corresponding to 
reversible binding reactions with the promoter of x,. Let c, =y,(Pi) 
denote the vector of concentrations of complexes with the 
promoter of x, at the quasi-steady state, obtained by setting 
0 = A i r i (p i ,c i ). 

We first define Hi(pi) as follows: 



y"=o 



(20) 



Table 1. Typical range of macroscopic parameters in Escherichia coli. 





Parameter 


Symbol 


Range 


Unit 


Reference 


Production rate constant 


n 


0-10000 


hr- 1 


[1], [49], [50] 


Dissociation constant 


k 


1-10 


nM 


[22], [52], [53] 


DNA concentration 


1 


1-500 


nM 


[51], [54], [56], [56] 


Protein decay rate 


d 


0.3-1.2 


hr" 1 


[1], [50] 



doi:l 0.1 371 /joumal.pcbi.1 003486.T.001 
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where yuiPi) is the j th entry in y ,-(/>,) and Q is the number of 
complexes with the promoter of x,-. 

Next, define the matrix *F; as follows: it has as many columns as 
the number of complexes formed with the promoter of x, , and as 
many rows as the number of parents of x, : 



C,,i C,-, 2 



Pi.l 

P/,2 



such that its (j,k) element is m if the j l parent of x, is bound as an 
ra-multimer in c,-^ (m = 0 if the j th parent is not bound). Finally, for 
nodes having parents, define the retroactivity Riipi) of node x, as 



dpi 



(21) 



For the most common binding types, Hj(pi) and Riipi) are given 
in Figure 3. For details on their derivation, see Supporting Text 
S3. 

Supporting Information 

Figure SI Mixing retroactivity can be avoided by 
introducing an extra node. Rectangles represent promoter 
regions, arrows denote coding regions. Promoters are regulated by 
TFs expressed from coding regions of the same color. (A) The 
production of is regulated by two TFs: x,- from the module, and 



Xj from the context. If the binding of x, is not completely 
independent from that of %j, the mixing retroactivity M of the 
context is non-zero. As a result, the dynamics of the context can 
suppress that of the module by (12). (B) One possible solution to 
obtain zero mixing retroactivity M is to introduce an extra input 
node x* in the context, so that parents from the module and from 
the context are not mixed. In particular, replace the promoter of 
Xfc with one that is regulated by x*. As a result, parents from the 
module and from the context are not mixed anymore, yielding 
M = 0, in the meantime, X/t still integrates the signal coming from 
x, (through x*) and from %j. 
(EPS) 

Text SI ODE model of the system in Figure 2 together 
with the parameter values used for simulation. 

(PDF) 

Text S2 Appendix containing the Theorems and Prop- 
ositions together with the corresponding proofs. Subsec- 
tions include the following: (1) Isolated module without input; (2) 
Isolated module with input; (3) Interconnection of modules; and (4) 
Bounding the difference between the trajectories of an isolated and 
a connected module. 
(PDF) 

Text S3 Derivation of Hi(p\) and Riipi) for the most 
common binding types presented in Figure 3. 

(PDF) 
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